*                           (-1)**n * e_n(r3) * r4**(2n+1)
*     sum for n=0 to inf of ------------------------------
*                                  2n+1 * exp(r3)
*     used for formula (15) and (16) of Litkouhi and Beck (1982);
*     e_n from Abramowitz and Stegun p262 eq6.5.11
      DOUBLE PRECISION FUNCTION RFLBSU(R3,R4)
      IMPLICIT NONE
      DOUBLE PRECISION R3, R4

      DOUBLE PRECISION DACC
      INTEGER          NORDMX
      PARAMETER (NORDMX=100)
      PARAMETER (DACC=1.E-6)
      DOUBLE PRECISION D3, D4, DSIGN, D3N, DENR3, D44, D42NP1, D2NP1,
     &                 DERM3, D0, DSUM, DI
      INTEGER          I

      IF (R3.GT.50.) THEN
         RFLBSU = 0.
      ELSE
         D3=DBLE(R3)
         D4=DBLE(R4)
         DSIGN  = 1.
         D3N    = 1.
         DENR3  = 1.
         D44    = D4**2
         D42NP1 = D4
         D2NP1  = 1.
         DERM3  = EXP(-D3)
         D0=D42NP1*DERM3
         DSUM=D0
         D0 = ABS(D0)*DACC
         DO I=1,NORDMX
            DSIGN = -DSIGN
            D3N   = D3N*D3
            DENR3 = DENR3+D3N/DBLE(I)
            D2NP1 = D2NP1+2.
            D42NP1= D42NP1*D44
            DI    = DSIGN*DENR3*D42NP1/D2NP1*DERM3
            DSUM  = DSUM+DI
            IF (ABS(DI).LE.D0) EXIT
         END DO
         RFLBSU=DSUM
      END IF

      RETURN
      END
